library(data.table)
library(Seurat)
library(SeuratDisk)
library(tidyverse)
library(tidyseurat)
source('00_util_scripts/mod_seurat.R')

all_marker <- c("IGKC","IGLC2","IGLC3","IGLC6","IGLC7","IGHA1","IGHA2","IGHG1","IGHG2","IGHG3","IGHG4", "IGHD","IGHE","IGHM","CD79A","CD79B","TRBC1","TRBC2","TRAC","TRDC","TRGC1","TRGC2","CD3D","CD3E","CD3G","CD247", "CD19","CD5")

# load smart-seq TPM data ---------
data <- fread('CRC-I/data/GSE146771_CRC.Leukocyte.Smart-seq2.TPM.txt.gz')

data <- column_to_rownames(data, 'V1')

II_smart <- c('P0104','P0411','P1212')
IT_smart <- c('P0305','P0720','P0825')

# load metadata
smrt_meta <- fread("CRC-I/data/GSE146771_CRC.Leukocyte.Smart-seq2.Metadata.txt.gz")

data[1:5,1:5]

## create seurat ------------
sobj <- data |>
  CreateSeuratObject(min.cells = 3,
                     min.features = 200,
                     names.field = 3)

# assign I232T genotype in smrt_meta
sobj <- sobj |>
  mutate(CellName = .cell,
         genotype = case_when(
           orig.ident %in% II_smart ~ 'II',
           orig.ident %in% IT_smart ~ 'IT',
           .default = 'NA'
         )) |>
  left_join(smrt_meta)

sobj$mitoRatio <- PercentageFeatureSet(object = sobj,
                                           pattern = "^MT-")

sobj |>
  VlnPlot('mitoRatio')

disease <- 'smart2-crc'

get_abundance_sc_wide(sobj, all_marker) %>%
  right_join(as_tibble(sobj)) %>%
  select(!matches("PC_|harmony")) %>%
  fwrite(str_glue("DE_cells/results/{disease}_de_meta.csv.gz"))

# save genotyped data
write_rds(sobj, 'DE_cells/data/CRC-smrt2-all.rds')

# load 10x TPM data ------------
data <- fread("CRC-I/data/GSE146771_CRC.Leukocyte.10x.TPM.txt.gz")

data[1:5, 1:5]

data <- column_to_rownames(data, 'V1')

II_10x <- c('P0410','P0323','P0408','P1026','P0104')
IT_10x <- c('P0123','P0202','P0613','P1025','P0305')

# the cell naming: P_N_P0410_01391, first letter is always 'P', second letter is N/P/T representing tissue source.

# load metadata
tenx_meta <- read_delim("CRC-I/data/GSE146771_CRC.Leukocyte.10x.Metadata.txt.gz")

sobj <- CreateSeuratObject(data,
                           names.field = 3,
                           min.cells = 3,
                           min.features = 200)

# assign I232T genotype in tenx_meta
sobj <- sobj |>
  mutate(CellName = .cell,
         genotype = case_when(
           orig.ident %in% II_10x ~ 'II',
           orig.ident %in% IT_10x ~ 'IT',
           .default = 'NA'
         )) |>
  left_join(tenx_meta)

sobj$mitoRatio <- PercentageFeatureSet(object = sobj, pattern = "^MT-")

sobj |>
  VlnPlot('mitoRatio')

disease <- '10x-crc'

# save 10x seurat file
get_abundance_sc_wide(sobj, all_marker) %>%
  right_join(as_tibble(sobj)) %>%
  select(!matches("PC_|harmony")) %>%
  fwrite(str_glue("DE_cells/results/{disease}_de_meta.csv.gz"))

# save genotyped data
write_rds(sobj, 'DE_cells/data/CRC-10x-all.rds')

# Chen 2021 data -------
# convert h5ad to h5seurat
list.files(path = 'DE_cells/data/cancer/', full.names = TRUE, pattern = 'h5ad') |>
  walk(Convert, dest = 'h5seurat')

Convert("DE_cells/data/cancer//VUMC_HTAN_VAL_EPI_V2.h5ad", dest = 'h5seurat')

## nonepi ---------
nonepi_sce <- zellkonverter::readH5AD('DE_cells/data/cancer/VUMC_HTAN_VAL_DIS_NONEPI_V2.h5ad',
                             reader = 'R',
                             verbose = T)

mymeta <- as_tibble(nonepi_sce)

nonepi <- nonepi_sce |>
  assay() |>
  CreateSeuratObject(min.cells = 3, min.features = 200) |>
  left_join(mymeta)

nonepi$mito_ratio <- nonepi |> PercentageFeatureSet('MT-')

# empirical cumulative distribution func
nonepi |>
  ggplot(aes(mito_ratio)) + stat_ecdf()

# 10.6k to 9.6k cells
nonepi <- filter(nonepi, mito_ratio < 40)

## discovery epi -----------
discov_epi <- LoadH5Seurat('DE_cells/data/cancer/VUMC_HTAN_DIS_EPI_V2.h5seurat', assays = 'counts')

discov_epi <- discov_epi |>
  UpdateSeuratObject() |>
  DietSeurat()

meta_discov <- as_tibble(discov_epi)

discov_epi <- GetAssayData(discov_epi) |>
  CreateSeuratObject(min.cells = 3, min.features = 200) |>
  left_join(meta_discov)

discov_epi$mito_ratio <- discov_epi |>
  PercentageFeatureSet(pattern = '^MT-')

discov_epi |>
  ggplot(aes(mito_ratio)) + stat_ecdf()

# 65k to 37k cells
discov_epi <- filter(discov_epi, mito_ratio < 40)

## validation epi ---------
valida_epi <- LoadH5Seurat('DE_cells/data/cancer/VUMC_HTAN_VAL_EPI_V2.h5seurat', assays = 'counts')

valida_epi <- valida_epi |>
  UpdateSeuratObject() |>
  DietSeurat()

meta_valida <- as_tibble(valida_epi)

valida_epi <- GetAssayData(valida_epi) |>
  CreateSeuratObject(min.cells = 3, min.features = 200) |>
  left_join(meta_valida)

valida_epi$mito_ratio <- valida_epi |>
  PercentageFeatureSet(pattern = '^MT-')

valida_epi |>
  ggplot(aes(mito_ratio)) + stat_ecdf()

# 57k to 50k cells
valida_epi <- filter(valida_epi, mito_ratio < 40)

## abnormal epi ---------
abnorm_epi <- LoadH5Seurat('DE_cells/data/cancer/VUMC_ABNORMALS_EPI_V2.h5seurat', assays = 'counts')

abnorm_epi <- abnorm_epi |>
  UpdateSeuratObject() |>
  DietSeurat()

meta_abnorm <- as_tibble(abnorm_epi)

abnorm_epi <- GetAssayData(abnorm_epi) |>
  CreateSeuratObject(min.cells = 3, min.features = 200) |>
  left_join(meta_abnorm)

abnorm_epi$mito_ratio <- abnorm_epi |>
  PercentageFeatureSet(pattern = '^MT-')

abnorm_epi |>
  ggplot(aes(mito_ratio)) + stat_ecdf()

# 23k to 18k cells
abnorm_epi <- filter(abnorm_epi, mito_ratio < 40)

## merge 4 subsets ----------
sobj <- c(discov_epi, valida_epi, abnorm_epi, nonepi) |>
  map2(c('discov_epi', 'valida_epi', 'abnorm_epi', 'nonepi'), AddMetaData, 'h5ad') |>
  purrr::reduce(merge)

sobj <- sobj |>
  tidyseurat::select(-c(Tumor_Type:last_col()))

# meta data
chen_meta <- read_csv('DE_cells/data/cancer/chen2021-meta.txt', col_names = FALSE) |>
  pull(X1) |>
  matrix(ncol = 14, byrow = TRUE) |>
  as_tibble() |>
  dplyr::select(1,3,14) |>
  set_names(c('specimen','sample','tissue'))

chen_meta |>
  count(tissue)

sobj <- sobj |>
  tidyseurat::left_join(chen_meta, by = c('HTAN.Specimen.ID' = 'specimen'))

sobj <- sobj |>
  filter(tissue != 'Not Otherwise Specified')

Idents(sobj) <- 'tissue'

VlnPlot(sobj, 'mito_ratio', pt.size = 0, group.by = 'tissue')

sobj <- sobj |>
  NormalizeData() |>
  FindVariableFeatures() |>
  ScaleData() |>
  RunPCA()

PCAPlot(sobj, group.by = 'tissue')

slice_sample(sobj, n = 100000) |> PCAPlot(group.by = 'tissue')

sobj <- harmony::RunHarmony(sobj, 'tissue')

sobj <- sobj |>
  RunUMAP(reduction = "harmony", dims = 1:20) |>
  FindNeighbors(reduction = "harmony", dims = 1:20) |>
  FindClusters()

sobj |> UMAPPlot(raster = F)

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj <- sobj |>
  mark_cell_type_singler(ref = hpca,
                         new_label = 'hpca_label')

sobj |> UMAPPlot(raster = F, group.by = 'hpca_label')

disease <- 'chen2021-crc'

write_rds(sobj, 'DE_cells/data/cancer/chen2021-mito40.rds')

sobj |>
  ggplot(aes(fill = hpca_label, x = tissue)) + geom_bar(position = 'fill')

# annotate DE cell features -----
get_abundance_sc_wide(sobj, all_marker) %>%
  right_join(as_tibble(sobj)) %>%
  dplyr::select(!matches("PC_|harmony")) %>%
  data.table::fwrite(str_glue("DE_cells/results/{disease}_de_meta.csv.gz"))

# AH-cOME ----------
disease <- 'AH_cOME'

get_abundance_sc_wide(sobj, all_marker, slot = 'count') %>%
  right_join(as_tibble(sobj)) %>%
  select(!matches("PC_|harmony")) %>%
  fwrite(str_glue("DE_cells/results/{disease}_de_meta.csv.gz"))
